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Abstract. We present calculations of the spatial and spectral distribution of the radio emission from a wide 
WR+OB colliding-wind binary system based on high-resolution hydrodynamical simulations and solutions to the 
radiative transfer equation. We account for both thermal and synchrotron radio emission, free-free absorption in 
both the unshocked stellar wind envelopes and the shocked gas, synchrotron self-absorption, and the Razin effect. 
To calculate the synchrotron emission several simplifying assumptions are adopted: the relativistic particle energy 
density is a simple fraction of the thermal particle energy density, in equipartition with the magnetic energy 
density, and a power-law in energy. We also assume that the magnetic field is tangled such that the resulting 
emission is isotropic. The applicability of these calculations to modelling radio images and spectra of colliding- 
wind systems is demonstrated with models of the radio emission from the wide WR+OB binary WR 147. Its 
synchrotron spectrum follows a power-law between 5 and 15 GHz but turns down to below this at lower and 
higher frequencies. We find that while free-free opacity from the circum-binary stellar winds can potentially 
account for the low-frequency turnover, models that also include a combination of synchrotron self-absorption 
and Razin effect are favoured. We argue that the high-frequency turn down is a consequence of inverse-Compton 
cooling. We present our resulting spectra and intensity distributions, along with simulated MERLIN observations 
of these intensity distributions. From these we argue that the inclination of the WR 147 system to the plane of 
the sky is low. We summarise by considering extensions of the current model that are important for models of 
the emission from closer colliding wind binaries, in particular the dramatically varying radio emission of WR 140. 
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1. Introduction 

Observations of early-type stars have revealed that they 
can be sources of both thermal and synchrotron radio 
emission. Thermal emission typically exhibits a bright- 
ness temperature ~ 10 4 K and a spectral index a ~ +0.6 
(S u oc v a ) at centimetre wavelengths, as expected from 
steady-state, radially symmetric winds (Wright & Barlow 
1975). In contrast, synchrotron emission is characterized 
by high brightness temperatures (> 10 6 K) and flat or neg- 
ative radio spectral indices. In addition to magnetic fields, 
synchrotron emission requires a population of relativistic 
electrons, which are widely thought to be accelerated in 
shocks. For single stars, shocks arise due to wind instabili- 
ties e.g. Lucy & White (1980), and propagate through the 
wind, while for massive binary systems stationary shocks 
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arise where the winds of the two stars collide e.g. Eichler 
& Usov (1993). Magnetic field compression within the 
shocked gas is another potential acceleration mechanism 
(Jardine et al. 1996). 

Spatially resolved observations of the WR+OB bi- 
nary systems WR 146 (Dougherty et al. 1996, 2000) and 
WR147 (Moran et al. 1989; Churchwell et al. 1992; 
Williams et al. 1997; Niemela et al. 1998) have dramat- 
ically confirmed the colliding-wind binary (CWB) model 
in these objects. In both WR 146 and WR 147 the thermal 
emission is coincident with the position of the WR star, 
while the synchrotron emission arises between the binary 
components at a position consistent with the pressure bal- 
ance of the two stellar winds. This model is supported fur- 
ther by the dramatic variations of the synchrotron radio 
emission in the 7.9-year WR+OB system WR140, which 
are clearly modulated by the binary orbit e.g. Williams 
et al. (1990a); White & Becker (1995). 
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Such results prompted van der Hucht et al. (1992) to 
suggest a CWB origin for all synchrotron emission from 
WR stars, for which there is now strong observational sup- 
port (Dougherty & Williams 2000). It is possible that this 
is also the case in O-type stars with synchrotron radio 
emission, e.g. Cyg OB2 #5 (Contreras et al. 1997), though 
there remain several apparently single stars that obviously 
do not fit a CWB interpretation. 

To date, modelling has been restricted to the radiom- 
etry from such systems. At frequency v, the observed 
flux (5° bs ) is related to the thermal flux (5** h ), the syn- 
chrotron emission arising from the wind-wind collision re- 
gion (Sl yn ), and the free- free opacity (r^) of the circum- 
system stellar wind envelope by 

S° hs = Sf + S7 n e-< (1) 

It is typically assumed that the stellar wind envelope is ra- 
dially symmetric and that the synchrotron emission arises 
from a point source at the stagnation point of the two stel- 
lar winds (see Williams et al. 1990a; Chapman et al. 1999; 
Skinner et al. 1999; Monnier et al. 2002, for examples). 
In this manner, relatively simple analytical solutions to 
the radiative transfer equation can be obtained. Although 
these models may successfully recover the radiometry, a 
single-valued free-free opacity determined along the line- 
of-sight to a point-like synchrotron emission region is an 
over-simplification, especially considering the extended re- 
gion of synchrotron emission observed from both WR 146 
and WR 147. Furthermore, in a CWB the assumption of 
radial symmetry fails in the collision zone. 

So far, no attempts have been made to construct syn- 
thetic radio images based on more realistic density and 
temperature distributions. This is surely imperative given 
the advent of spatially resolved observations where direct 
comparison between models and observations may be ex- 
pected to increase our understanding of this phenomenon. 
In making the first steps toward addressing this situation, 
we have calculated the free-free and synchrotron radio 
emission arising from an early-type binary system under 
various simplifying assumptions. Our approach extends 
and improves on previous work by including the ability to 
simulate both the free-free and synchrotron emission and 
absorption from the stellar winds and wind-wind collision 
region based on 2D, axis-symmetric hydrodynamical sim- 
ulations of the density and pressure distribution. Various 
synchrotron cooling mechanisms can also be incorporated 
as required. An arbitrary inclination angle for the line of 
sight can be specified, and the radiative transfer equation 
is then solved to generate synthetic images and spectra. 

The layout of the remainder of this paper is as follows. 
In Sec. 2 we describe our model in more detail and in Sec. 3 
we perform a parameter study using a "standard" model 
of a very wide CWB. We examine the influence of syn- 
chrotron self-absorption and the Razin effect, which have 
not hitherto received much attention, along with system 
inclination and binary separation. The application of our 
model to observations of WR 147 is described in Sec. 4. 



In Sec. 5 we summarize and note future directions. In 
Appendix A, the geometry of the ray-tracing technique 
for solving the radiative transfer equation is described. 

2. Modelling the radio emission from CWBs 

To calculate the radio emission from an early-type bi- 
nary, we must first compute the structure of the wind- 
wind collision. This is most readily achieved through the 
use of a hydrodynamical code, and we use VH-1 (Blondin 
et al. 1990), a Lagrangian remap version of the piecewise- 
parabolic method e.g. Pittard & Stevens (1997). In this 
paper we are concerned particularly with wide systems, so 
it is assumed that the spherically symmetric stellar winds 
are accelerated instantaneously to terminal speeds. This 
is reasonable given that the winds attain terminal veloc- 
ity at a small percentage of the radial distance from the 
stars to the wind collision region. Distortion of the colli- 
sion zone by orbital motion, most especially close to the 
shock apex, is negligible so axis-symmetry is also assumed. 
We use WR or solar abundances for the stellar winds, as 
appropriate. Such models have been successfully applied 
to the analyses of X-ray observations from CWBs (see 
Pittard & Corcoran 2002; Pittard et al. 2002, and refer- 
ences therein). 

Once a suitable hydrodynamic solution has been ob- 
tained, the necessary information is read into our radiative 
transfer ray-tracing code, and appropriate emission and 
absorption coefficients for each cell on the 2D grid are de- 
termined. A synthetic image on the plane of the sky is then 
generated by solving the radiative transfer equation along 
suitable lines of sight through the grid (sec Appendix A 
for more details). The following sections detail the calcu- 
lation of the emission and absorption properties for each 
hydrodynamic cell. 

2.1. Thermal emission and absorption 

The thermal emission (e^) and absorption (a^) coeffi- 
cients at a frequency v are given by Rybicki & Lightman 
(1979) as 

el = 6.8 x lO-^Z^^T-^e-^/^gs, (2) 
oF v = 3.7 x ltfT- l / 2 Z 2 n e n lV -*{l - e- hv / kT )g s , (3) 

where Z is the ionic charge, n e and rii are electron and 
ion number densities, T is the temperature of the gas, and 
gs is a velocity averaged Gaunt factor (Hummer 1988). 
The densities n e and rij are determined from the cell den- 
sity, composition and ionization, where the ionization is 
specified £ls Si function of cell temperature. The appro- 
priate ionization of the elemental species for a particular 
temperature regime is determined and the absorption and 
emission coefficients are then evaluated for several ionic 
charges. 

The accuracy of the "thermal" code was verified by 
using a simple spherically symmetric, isothermal stellar 
wind model with an r~ 2 radial density distribution. The 
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resulting fluxes and free-free opacities matched those de- 
rived from the analytical expressions given in Wright & 
Barlow (1975). 

2.2. Synchrotron emission 

When non-relativistic particles are accelerated in a mag- 
netic field, emission occurs at frequencies corresponding 
to harmonics of the gyration frequency of the particle. 
This is known as cyclotron emission, and most of the radi- 
ated power occurs at the gyration frequency. On the other 
hand, if the particles are relativistic, the emission blurs 
into a continuum rather than a series of delta functions 
and emission may occur at a frequency many orders of 
magnitude greater than the gyration frequency. The syn- 
chrotron emission from a single electron is (cf. Rybicki & 
Lightman 1979) 



P(u) 



V3q 3 B 



sin a 



F{v/u c ) 



(4) 



where q is the particle charge, B is the magnetic field 
strength, a is the pitch angle of the particle relative to the 
direction of the B-field, m is the mass of the particle, and 
F{y/v c ) is a dimensionless function describing the total 
power spectrum of the synchrotron emission, with v c the 
frequency where the spectrum cuts off, given by 



3-f 2 qB sin a 
Aitmc 



(5) 



Values for F(v/v c ) are tabulated in Ginzburg & 
Syrovatskii (1965). Throughout the remainder of this pa- 
per we assume that sin a = 1. The emissivity per unit 
volume from a mono-energetic distribution of electrons is 
then simply the product of Eqn. 4 and the number density 
of the radiating particles. 

In our model we suppose that the relativistic elec- 
trons are created by lst-order Fermi acceleration at the 
shocks where the winds collide. For a strong shock with 
a compression ratio of 4, test particle theory predicts 
that accelerated electrons will have an energy distribu- 
tion with a power-law index of 2 (see references in Eichlcr 
& Usov 1993). Hence, the number of relativistic electrons 
per unit volume with energies between 7 and 7 + dj is 
N{^)d^ — C"f~ p d"/, where 7 is the Lorentz factor and 
p = 2. More recent nonlinear calculations yield a simi- 
lar energy spectrum for a wide range of shock parameters 
e.g. Ellison & Eichler (1985). A useful discussion of shock 
acceleration can be found in Ellison & Reynolds (1991). 

For such a power-law distribution of electrons, the 
total synchrotron emission at frequency v is given by 
(Rybicki & Lightman 1979) 



V3q 3 CB 
yV) mc 2 (p+l) 



p 19 

- H 

4 12 



£__L 

4 12 
\ 3qB ) 



(6) 



where T is the Gamma function. Eqn. 6 is only strictly 
valid in the synchrotron limit, when the frequency of emis- 
sion is between approximately a few times the gyration 
frequency (vb) and v c (defined in Eqn. 10), but we re- 
strict this investigation to frequencies much greater than 
vb e.g. for our model of WR 147 vb < 10 4 Hz at the shock 
apex. 

To evaluate Eqn. 6, we must determine the appropriate 
value for the normalization constant C, and the magnetic 
field strength B. Since our hydrodynamical simulations do 
not provide direct information on the magnetic field or rel- 
ativistic particle distribution, we follow the standard pro- 
cedure in such cases (e.g. Chevalier 1982; Mioduszewski 
et al. 2001): that is, to assume that the magnetic energy 
density Ub, and the relativistic particle energy density 
U IC \, are proportional to the thermal particle internal en- 
ergy density C/ t h- Then 



U B = B 2 /8k = Cst/th, 

U le \ = J n(^)-fmc 2 d-f = Crel^th, 



(7) 
(8) 



where (b and Crci are constants which are typically set 
to around 1% (Mioduszewski et al. 2001), and {/ t h = f-P 
where P is the gas pressure and the adiabatic index of 5/3 
is assumed, as for an ideal gas. The common assumption 
of equipartition corresponds to Qb — Col = C 

As noted in Mioduszewski et al. (2001), there is some 
justification for such a relation between the energy den- 
sities: the thermal energy density is higher in the post- 
shock than in the pre-shock gas, and this is also where the 
magnetic field is likely to be enhanced, and the particles 
accelerated to relativistic energies. However, such argu- 
ments ultimately need to be replaced by detailed physical 
processes. It is well established that the shock acceleration 
of protons is very efficient, but how much energy is trans- 
ferred to electrons in shocks is much less certain. We note 
that the relativistic electron energy is unlikely to be more 
than 5% of the total available shock energy, and could con- 
ceivably be much less (Blandford & Eichler 1987; Ellison 
& Reynolds 1991; Eichler & Usov 1993). 

For p = 2, 



C 



rcl 



m e c 2 lir/n 



(9) 



where we have assumed 1 < 7 < 7 max and 7 max specifies 
the maximum energy of the relativistic electrons. Since 
it is assumed in Eqn. 6 that 7 max — > 00, the synchrotron 
emission extends to infinitely high frequencies with a spec- 
trum cx v~ 5 . In reality the synchrotron flux will depart 
from this relationship at frequencies near v Cl where v c is 
the frequency of peak emission from a particle of energy 
E = 7 max TOC 2 , given by (Rybicki & Lightman 1979) 



3 il^qB 

4n 



mc 



(10) 



In CWBs we expect 7 max to be set by the balance between 
energy gain by Fermi acceleration and various cooling 



4 



S.M. Dougherty et al.: Radio emission models of Colliding- Wind Binary Systems 



processes, which include inverse- Compton (hereafter IC) 
cooling, synchrotron decay, and ion-neutral wave damp- 
ing. Ion-neutral damping does not play a dominant role in 
CWBs (Eichler & Usov 1993), and is not discussed further 
here. The relevant processes and timescales are discussed 
in the following subsections. 

2.3. Synchrotron self-absorption 

We include synchrotron self-absorption (hereafter SSA) 
in our models. For a power-law distribution of particle 
energies, the absorption coefficient for SSA is (Rybicki & 
Lightman 1979) 



8irm \ 2irm 3 c 5 



3q 



p/2 



CB ip+2)/2 



r(^)r(>±^),^. (ii) 

where the symbols have their previously defined meanings. 



2.4. The Razin effect 

When relativistic charges are surrounded by a plasma 
(as opposed to existing in a vacuum), the beaming effect 
that characterizes synchrotron radiation is suppressed. In 
essence, the refractive index of the medium reduces the 
Lorentz factor of the charge to 



7 



v/1 + 7^0 2 /^ 2 



(12) 



where the plasma frequency vq — U q 2 n/irm. Hence syn- 
chrotron emission by relativistic charges in the medium 
is possible only when 7' >> 1, and a low frequency 
cut-off occurs. Quantitative details of this process, which 
is known as the Tsytovitch-Eidman-Razin effect, can be 
found in Hornby & Williams (1966). The characteristic 
cut-off frequency is given by 



vr = 20-^. 



(13) 



Since the emission decreases exponentially at frequencies 
low compared with vr, and there is a noticeable effect 
at frequencies greater than j/r, e.g. there is a 10% reduc- 
tion in flux at v = IQisr, we approximate this reduction 
by multiplying the synchrotron flux from our model by a 
factor e~ VR l v . 



2.5. The low and high-energy cutoff in the relativistic 
particle distribution 

The nature of the synchrotron spectrum is determined 
by the underlying relativistic particle energy distribution. 
The low and high energy cutoffs in the energy distribu- 
tion are particularly important. At low energies, the rel- 
ativistic particle distribution can be potentially affected 
by coulombic collisions with thermal ions. However, at the 



densities and temperatures in the wind-collision regions of 
the wide systems considered in this paper this process is 
insignificant. The high energy cutoff is determined at the 
point where the rate of energy gain balances the rate of 
energy loss. Radiation losses due to synchrotron emission 
and IC scattering are in the same ratio as the magnetic 
field energy density to the photon energy density, 



Psvnc. Ub 



sync 
^compt 



ph 



(14) 



As noted in Sec. 2.2, it is expected that (b is no more than 
0.05, which suggests that P syn c < -Pcompt m the wide CWB 
systems under consideration, and IC cooling dominates 
over synchrotron losses. 

Electrons are rapidly accelerated to relativistic ener- 
gies at the shocks bounding the wind-collision zone, and 
attain an energy spectrum specified by p = 2 and 7 ma x 
(where the rate of energy loss by IC cooling balances the 
rate of energy gain from 1st order Fermi acceleration). 
They are then advected out of the shock and (for 7 < 10 5 ) 
are effectively frozen into the post-shock flow (White 
1985). If the characteristic flow time out of the system 
is shorter than the timescale for IC cooling, the dominant 
cooling process is adiabatic expansion. This is the case 
for wide CWBs and low 7, and since adiabatic cooling is 
treated by the hydrodynamic code, we assume that (3 and 
Croi are spatially invariant in the post-shock flow. For high 
7, IC losses can be very rapid even in the widest systems, 
as we show is the case for WR147 in Sec. 4. IC cooling 
will be an even more important consideration in closer 
CWBs e.g. WR140 and shorter period systems, where it 
could prevent the acceleration of electrons to even fairly 
modest values of 7. In these first attempts at modelling 
the radio emission from wide CWBs we have not included 
IC cooling explicitly in the code, but where relevant we 
discuss its impact on the resulting spectra e.g. Sees. 3.3 
and 4. This will be treated more fully in Pittard et al. (in 
preparation) . 

2.6. Some simplifying assumptions 

For simplicity in this first investigation, we assume that 
the post-shock ion and electron temperatures equalize 
very rapidly. While the equilibration timescale can be sig- 
nificant compared with the flow timescale e.g. WR 140 - 
see Usov (1992); Zhekov k Skinner (2000), the exact sit- 
uation is currently unclear with the importance of pos- 
sible electron heating mechanisms in collision-less shocks 
still being debated. We anticipate that our simplification 
will have little bearing on the resulting radio emission be- 
cause the synchrotron emission from the wind-collision 
zone dominates the thermal free-free emission. While it 
may affect the efficiency with which the accelerated elec- 
trons are extracted from the thermal pool, this conjecture 
is impossible to test currently. Instead its effect is simply 
folded into the value of £. Also, we assume that the gas 
is in ionization equilibrium, which is not most likely the 
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Fig. 1. Spectra from our standard model with 0° inclination 
- free-free flux (dashed), synchrotron flux (dotted), and the 
total flux (solid) . The effect of free-free absorption on the syn- 
chrotron emission in this wide system is only just evident below 
300 MHz. The spectral index of the thermal spectrum is +0.6, 
and the optically thin part of the synchrotron spectrum is —0.5. 
Neither the Razin effect nor SSA are treated in this model. 

case immediately post-shock, especially in wide systems. 
However, we believe that this will not impact the results 
presented here since the bulk of the emission from the 
shocked gas is synchrotron, not free-free emission. 

One final assumption in our models is that the 
magnetic field in the post-shock gas is highly tangled. 
Conveniently, this allows us to treat the synchrotron emis- 
sion and absorption as isotropic, and means that it is pos- 
sible to use isotropic formulae to calculate the IC losses 
(Pittard et al. , in preparation). While the magnetic held 
from each star should be structured on large scales, and 
toroidal in shape at the large distances considered in the 
models presented here (see Eichler & Usov 1993, and ref- 
erences therein), one might expect finer-scale structure in 
the post-shock flow. However, such details are beyond the 
scope of this paper. 

3. Parameter Study 

3.1. A standard model 

Our initial investigations have been based around a stan- 
dard CWB model with the following parameters: M W r = 
2 x 1(T 5 Mq yr- 1 , M OB = 2 x 1(T 6 M Q yr" 1 , Uo0;WR = 
Woo.ob = 2000 kms" 1 , D = 2 x 10 15 cm, and T wind for 
both stars of 20 kK. These wind values are typical of WR 
and OB stars. At a frequency v — 5 GHz, the adopted 
binary separation D, is approximately lOx the radius of 
the Tff = 1 surface of the WR wind. For inclination an- 
gles ^ 0°, the lines of sight to the wind-collision zone are 
then optically thin, permitting investigation of the result- 
ing emission from the wind-collision region in the absence 
of strong free-free absorption from the circum-binary stel- 
lar wind envelope. We also assume solar abundances for 
the OB star and mass fractions X = 0, Y = 0.75, Z = 0.25 



appropriate for a WC-type star, and an ionization struc- 
ture of H+, Hc+ and CN0 2+ . We adopt C = 10~ 4 as our 
"canonical" value of U-B./Uth and place our model system 
at a distance of 1.0 kpc. Maximum values of temperature, 
density, and magnetic field occur at the shock apex and 
for these input parameters we find T max = 1.75 x 10 8 K, 
Pmax = 9x 10~ 19 gcuT 3 (n max = 4 x 10 5 cur 3 ), and 
B max = 6 mG. Though B is an equipartition value, fields 
of this order at the location of the shock extrapolate to 
surface fields on the OB star of <~ 30-300 G, using the 
magnetic field structure given in Eichler & Usov (1993) 
and taking the typical rotation speed of O-type stars to 
be ~ O.Iwoo^ob (Conti & Ebbets 1977). Such surface fields 
are consistent with the few observations of the fields in 
these types of stars (Donati et al. 2001, 2002). For all of 
the models presented in this paper we fix 7 max — 4000. 
Since we are using Eqn. 6 to calculate the synchrotron 
luminosity, fixing 7 max only has a weak impact on the 
constant C <~ 1/ln 7 max . 

The radiomctry and intensity distributions obtained 
from our standard model are shown in Figs. 1 and 2. In 
these first examples we only consider the impact of free- 
free absorption: SSA and the Razin effect arc examined in 
the next section. Around 1 GHz the synchrotron emission 
from the wind-wind collision is comparable to the total 
free-free emission, while at 22 GHz the emission is largely 
from the two stellar winds. Below 1 GHz, the synchrotron 
emission dominates the total flux and the beginnings of a 
turnover seen at ~ 300 MHz in the synchrotron spectrum 
is due to free-free absorption from the unshocked stellar 
winds. Above 300 MHz, the synchrotron spectrum is a 
power-law with a spectral index equal to (1 — p)/2 = —0.5 
for p = 2. The thermal spectrum is also a power- law with 
a spectral index of +0.6, as expected for a fully ionised, 
isothermal stellar wind envelope with a r~ 2 radial density 
distribution. 

Although we have deliberately avoided modelling a 
specific system in this first analysis, these initial results 
already bear similarities to observations e.g. WR 147 - see 
Sec. 4. The total flux, the spectral shape, and the spatial 
intensity distribution are of the correct order of magnitude 
and morphology. 

3.2. Razin effect and Synchrotron Self-Absorption 

In Figs. 3 and 4 the influence of the Razin effect and SSA 
on the synchrotron spectrum of our standard model is 
shown compared with the case where neither mechanism 
is active. Note that in the models shown, free-free opacity 
is negligible and has little influence on the spectra. It is 
clear that both the Razin effect and SSA can be impor- 
tant in quenching the low-frequency spectrum, with the 
relative impact of these two mechanisms in these models 
being strongly dependent on the value of (. For £ = 10~ 4 , 
the Razin effect is the dominant mechanism whereas for 
C = 10~ 3 , SSA dominates. From Eqs. 11 and 13, u v cx ( 2 
and v R oc (for fixed n e ) so we expect SSA to be in- 
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Fig. 2. Intensity distributions of our standard model at 0° in- 
clination at 1.6 GHz (top) and 22 GHz (bottom). Both im- 
ages have the same intensity scale, highlighting the relative 
importance of the emission from the wind-wind collision re- 
gion and the stellar winds as frequency varies (see also Fig. 17). 
Though not shown, the emission extends to the panel bound- 
aries. Neither the Razin effect nor SSA are included in this 
calculation. 



Frequency (GHz) 

Fig. 3. The influence of the Razin effect on synthetic syn- 
chrotron spectra using the standard model at 0° inclination 
for various values of £. The effect of free-free absorption by the 
circum-binary stellar wind envelope is negligible, just visible 
at around 200-300 MHz. No Razin (solid), C = 10~ 3 (dotted), 
10~ 4 (dashed), 10~ 5 (dash-dotted). All spectra are normalized 
to the 80 GHz synchrotron flux when no Razin effect is in- 
cluded. 
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Frequency (GHz) 

Fig. 4. The effect of SSA on synthetic synchrotron spectra us- 
ing the standard model at 0° inclination for various values of 
£. The effect of free-free absorption by the circum-binary stel- 
lar wind envelope is negligible, just visible at around 200-300 
MHz. No SSA (solid), C = 10" 4 (dotted), 10~ 3 (dashed), 10" 2 
(dash-dotted). All spectra are normalized to the 80 GHz syn- 
chrotron flux when no SSA is included. 



creasingly important as ( increases, with the Razin effect 
becoming more important for decreasing £. For our stan- 
dard model and reasonable values of 7 max , we find that the 
Razin effect dominates the low-frequency absorption when 
£ < 10 -4 , the Razin effect and SSA are roughly compara- 
ble when C ~ 10~ 3 , and SSA dominates when ( > 10~ 2 . 
Of course, we expect these values to be different for other 
model systems. 

For a given value of £, a different value of 7 max would 
affect the value of C (Eqn. 9) and hence the intrinsic syn- 
chrotron flux, P(v) (Eqn. 6), although the synchrotron 
spectrum would still extend to infinite frequency as a re- 



sult of adopting Eqn. 6. To obtain a certain intrinsic flux, 
only the combination C 7 ^ 4 / m 7max is required to be a spe- 
cific value and £ and 7 max can vary within this constraint. 
However, since a v (Eqn. 11) and i/r (Eqn. 13) depend 
on different powers of £ and 7 max , this degeneracy breaks 
down when SSA and the Razin effect are important. This 
means that the application of our model is dependent on 
the particular value of 7 max assumed. Nevertheless, since 
ln7max will vary by less than a factor of 4 for a reason- 
able range of 7 max , the strength of these absorption terms 
will not change a great deal. In fact, as a v cx CV m 7m ax 




1 10 100 

Frequency (GHz) 

Fig. 5. The effect of separation on synchrotron spectra using 
the standard model at 0° inclination and £ = 10~ 4 , for D — 
2, 5, 20, and 80 x 10 14 cm. Only the effect of free-free opacity 
is included. 

i.e. almost the same dependence as the intrinsic flux, its 
variation will be very small. As vr oc C~ 1//2 i the relative 
strength of the Razin effect is slightly more sensitive to 
our particular choice of 7 max , but not hugely so. Hence 
the curves in Figs. 3 and 4 are not strictly functions of 
only C, as there is a slight implicit dependence on 7 max as 
well. 

Figs. 3 and 4 also demonstrate that the spectral shape 
belies the underlying mechanism; the SSA spectra are 
power-law, whereas the Razin dominated spectra decay 
exponentially. Optically thick synchrotron emission actu- 
ally has a power-law spectrum with a slope of +2.5, but 
we find that the spectral slope in our calculations is ~ +1. 
This difference reflects the fact that the synchrotron emis- 
sion in our models is always a mixture of optically thick 
and thin emission. When the synchrotron emission from 
the shock apex is optically thick, a substantial amount 
of optically thin emission from the downstream flow con- 
tributes to the total non-thermal emission. 

3.3. Effect of Binary Separation 

The dramatic variations in the radio flux of the highly ec- 
centric WR+OB binary WR 140 point to the significance 
of binary separation to the observed radio emission. In 
WR140, the separation varies from <~ 2 AU at perias- 
tron to ~ 28 AU at apastron, a separation range over 
which different physical mechanisms determine the result- 
ing synchrotron spectrum. In this section, we perform cal- 
culations with varying binary separation to investigate the 
effect on the intrinsic synchrotron luminosity, and also how 
separation affects the relative importance of free- free opac- 
ity, the Razin effect, and SSA on the spectrum. We also 
discuss the impact of IC cooling as a function of separa- 
tion. 

In Fig. 5 we show the effect of free-free opacity on 
the spectra resulting from varying the separation between 



1 10 100 

Binary Separation (10 14 cm) 

Fig. 6. The effect of separation on the synchrotron luminosity 
for a power-law electron energy distribution in the optically 
thin regime. Shown are the data points at four frequencies: 1.6 
(solid), 5.0 (dotted), 8.3 (dashed), 14.6 GHz (dot-dashed). The 
slope of this log-log plot is -1/2, as expected. 
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Fig. 7. The effect of separation on the free-free opacity within 
the wind-wind collision region, using the standard model at 0° 
inclination and C = 10~ 4 , for D = 2, 5, and 20 x 10 14 cm. The 
intrinsic synchrotron emission is shown (solid) along with the 
spectra of the synchrotron emission emerging from the collision 
zone (dotted). 



2 x 10 cm and 80 x 10 cm, while the other parameters of 
our standard model were unchanged. The radius of optical 
depth unity for the WR wind is w 2 x 10 14 cm at 5 GHz. 
The increasing importance of free-free absorption as the 
separation decreases is clearly observed at the lower fre- 
quencies. As the stars move closer together the wind- wind- 
collision region moves closer to both the WR star and the 
companion, and is surrounded by increasingly dense gas 
which increases the line-of-sight opacity to the shock apex. 
For a stellar wind with an r~ 2 radial density distribution 
it can be easily shown that the line-of-sight opacity r, at 
frequency v through the wind is r cx £ _3 i/ -21 , where £ 
is the impact parameter. Since £ <x D for a given inclina- 
tion, the turnover frequency for a constant opacity value 
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Fig. 8. The effect of separation on the synchrotron emission 
using the standard model at 0° inclination and £ = 1CP 4 , for 
D = 20 x f0 14 (labelled "20"), 5 x f0 14 cm ("5") and D = 
2 x 10 14 ("2"): the intrinsic synchrotron emission i.e. no Razin 
or SSA (solid); only Razin effect included (dashed); only SSA 
included (dotted). The intrinsic synchrotron emission is shown 
to highlight the effect of separation on the synchrotron emission 
i.e. in the absence of the strong free-free absorption, as seen in 
Fig. 5. 

is v oc D~ 10//7 , in excellent agreement with the spectra in 
Fig. 5. 

With the parameters for our standard model, Fig. 5 
also reveals that the synchrotron luminosity increases as 
the separation decreases, due to the increased thermal en- 
ergy density in the collision region. From Eqn. 6, it can be 
deduced that the intrinsic synchrotron emission per unit 
volume P{v) oc (^ 3 / 4 n 3//4 ^ _1 / 2 for an electron power-law 
index p = 2, in the absence of SSA or the Razin effect. 
Since the post-shock density cx D~ 2 and the volume of the 
emitting region scales as D , the total synchrotron emis- 
sion from the entire wind collision volume oc D~ x l 2 v~ x l 2 . 
This is illustrated in Fig. 6. For comparison, we note that 
the total intrinsic X-ray emission from the wind-wind col- 
lision scales as D~ x in the optically thin, adiabatic limit 
(Stevens et al. 1992). 

In Fig. 7 the effect of free-free opacity from within 
the shocked gas is shown. As the binary separation is re- 
duced the free-free opacity in the shocked region becomes 
increasingly important, as density increases. Note that 
shocked gas temperature is independent of separation. As 
before, the free- free spectral turnover oc Z)~ 10 / 7 , in close 
agreement with the spectra shown. In this case, in systems 
such as WR 140 with periaston separation ~3x 10 13 cm, 
the free-free opacity of the shocked plasma will have an 
impact of the spectrum below ~ 10 GHz. 

The impact of changing separation on the importance 
of the Razin effect and SSA is shown in Fig. 8. The in- 
trinsic synchrotron spectrum is shown in order that the 
dramatic effects of the free-free opacity of both the un- 
shocked circum-binary stellar winds and the shocked gas, 
clearly seen in Fig. 5 and 7, are removed. As expected, 
the attenuation from both processes is more pronounced 




-100 -50 " 50 



Fig. 9. Intensity distributions at 1.6 GHz and 0° inclination for 
an orbital separation of 20 (top), 5 (middle), and 2 x 10 14 cm 
(bottom), corresponding to angular separations of 133, 33 and 
13.3 mas respectively. The crosses denote the positions of the 
stars, with the WR star located at (0,0). Each image has the 
same intensity scale and contours. The contours were chosen to 
highlight the relative brightness of the stellar wind of the WR 
star and the emission from the collision region. Though not 
shown, the emission extends to the boundaries of the panels 
(see also Fig. 17). 

as the separation decreases since the shock apex occurs 
in a higher density part of the stellar wind envelope. The 
Razin turnover frequency i>r oc nj B, and as B oc (for 
fixed Q, v R increases as n'. Because n oc D~ 2 , we find 
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that v R cx Z) . With D = 2 x 10 5 cm, our standard 
model gives vr = 1.3 GHz for emission at the shock apex. 
Away from the shock apex vr will be smaller. For our 
model with D = 5 x 10 14 cm, vr = 5 GHz. Both of these 
values are in excellent agreement with the frequency where 
the Razin effect causes a 1/e reduction in flux, as seen in 
Fig. 5. The opacity due to SSA is tssa ~ ct v roB, where 
tob is the distance between the shock apex and the OB 
star. Since a u cx D~ A v~^ and tob D (see Eqn. 15), this 
leads to fssA oc D^ 1 , which is in excellent agreement with 
the spectra shown in Fig. 5, where the SSA turnovers oc- 
cur at approximately 0.9, 3.5, and 9.0 GHz for separations 
of 2, 5 and 20 x 10 14 cm respectively. 

Also of interest is how the distribution of the syn- 
chrotron emission is affected by changing separation. This 
is shown in Fig. 9 for three different separations. In the 
middle and bottom images, the maximum intensity ac- 
tually occurs to either side of the axis of symmetry, this 
offset widening as the separation decreases. This is consis- 
tent with significant absorption along the line of sight to 
the shock apex, and slightly reduced absorption to the po- 
sitions of peak emission. In addition, a given flux intensity 
from the wind collision is apparent out to greater off-axis 
distances as the separation decreases. This is due to a com- 
bination of increased energy density in the collision zone 
and the importance of absorption from the surrounding 
winds. 

As mentioned in Sec. 2.5, IC scattering is an important 
cooling mechanism for relativistic electrons. To estimate 
when IC cooling needs to be considered, we determine the 
minimum value of 7 at which IC cooling can have an effect 
before the relativistic electrons adiabatically flow out of 
the system. The distance from the OB companion to the 
shock apex is given by 



roB = 



1 + 7?5 



-D. 



(15) 



Here, r) is the momentum ratio of the two stellar winds 
in the system and D is the binary separation. For our 
standard model, 77 = 0.1 and tob/D = 0.24. Assuming a 
strong shock, and that the emission is concentrated near 
the shock apex within an off-axis distance of TrroB/2 1 , the 
characteristic flow timescale for gas is ifl ow = 2irroB/v, 
where the post-shock velocity is v/A, with v being the 
relevant pre-shock velocity of either the WR or OB stellar 
wind. 



1 This is true in the wider systems, but in close systems 
strong absorption along the line of sight to the shock apex can 
create flux maxima to either side of the shock apex, as shown 
in Fig. 9. However, in such cases a substantial proportion of the 
mass flux (relativistic electron flux) through these regions will 
have passed through (originated in) the shock at some non- 
negligible off-axis distance. Therefore, the following argument 
should be good to first order, even in close systems. 



For relativistic particles with v/c w 1, the rate of en- 
ergy loss through IC cooling is (Rybicki & Lightman 1979) 



dE 
~dt 



(7t7 



IC 



3tt 



(16) 



where is the stellar luminosity and r is the distance 
of the shock from the star. Since the relative momentum 
fluxes of the WR and OB stellar winds normally cause 
the collision zone to be closer to the OB star, which is 
generally also the more luminous star in the binary, L* 
and r should be evaluated for the OB star. From Eqn. 16, 
the timescale for 50% energy loss through IC scattering is 



tic = 



37rm e c 2 rQ B 



(17) 



where we have assumed that to first order the average dis- 
tance of the relativistic electrons from the source of UV 
photons during this time is tob- If the distance from OB 
star (the dominant source of UV photons) to the relativis- 
tic electrons is such that an electron of energy jm e c 2 loses 
50% of its energy through IC cooling during the flow time, 
then t flow = ti c when 



_ 2a T L < .7fl ow 
3m e c 2 v 

or alternatively, when 



7flc 



3m e c 2 vr B 

2<7t£* 



(18) 



(19) 



If we assume further that the luminosity of the OB star 
is L* ~ 10 5 L Q , then from Eqn. 19 we find 7 w 500 for 
our standard model, which from Eqn. 10 gives a charac- 
teristic frequency v c w 6 GHz for B = 6 mG. Hence, for 
our standard model, the above order of magnitude esti- 
mate suggests that IC cooling starts to have an effect on 
the relativistic electrons with 7 > 7fl ow - In practice, this 
is likely to represent a lower limit. Firstly, it has been de- 
rived for electrons accelerated at the symmetry axis and 
which then flow to a distance of 7rroB/2 from the sym- 
metry axis. Electrons accelerated somewhat off-axis will, 
of course, take less time than ta ow to reach this off-axis 
distance. Moreover, the average distance to the OB star 
will be greater than ros- Therefore, we expect that IC 
cooling has only a significant effect for 7 greater than a 
few times 7fl ow - For our standard model, electrons with 
7 ~ 37fiow result in emission with a characteristic fre- 
quency of ~ 50 GHz. 

This estimate suggests that IC cooling is an important 
consideration in all systems, regardless of binary separa- 
tion, for sufficiently high frequencies (or 7). However, it 
is clearly most significant in shorter period systems, since 
7fi ow cx D and v c cx j 2 B cx D (assuming B cx 1/-D). 
Models with IC cooling included explicitly will be investi- 
gated in Pittard et al. (in preparation). 
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Fig. 10. The effect of free-free opacity as a function of inclina- 
tion angle on the synchrotron spectra of the standard model 
with = 1CP 4 . The inclination angle of the axis of symme- 
try to the plane of the sky for the various spectra is shown. 
Dashed lines indicate i < 0° (i.e. WR star in front), the solid 
line i = 0° (i.e. quadrature), and dotted lines i > 0° (i.e. OB 
star in front). 



3.4. Effect of inclination angle 

We now examine the effect on the observable synchrotron 
spectrum of changing the inclination of the axis of sym- 
metry (the line of centres) to the line of sight using our 
standard model. This is shown in Fig. 10, where only the 
effect of the free-free opacity of the circum-binary stel- 
lar wind envelope is being considered. Since our model is 
axis-symmetric, changing only the inclination covers any 
orientation of the system relative to the observer. The 
largest effects are seen at the lowest frequencies, where 
the free-free opacity is highest for a given path length 
through the circum-binary envelope (rg oc v~ 2A ). At an 
inclination angle i = 0°, the lines-of-sight are perpendic- 
ular to the axis of symmetry of the model; at i = —90° 
they pass first through the WR-star wind; and at i = +90° 
they pass first through the OB-star wind. The asymptotic 
half-opening angle of the WR-shock cone is 72°, so when 
i > 18°, lines of sight into the system can impinge the 
shocked region without first passing through unshocked 
circum-binary gas. 

The effects of changing free-free opacity with inclina- 
tion are clear. The size of the optically thick region of the 
circumbinary nebula goes as f~ - 7 for an r~ 2 radial den- 
sity gradient, and so as frequency increases the opacity 
is reduced along any given line-of-sight. In our standard 
model, it is evident that the cirumbinary nebula is opti- 
cally thin above a few GHz. The largest free-free attenua- 
tion occurs when the line of sight to the wind-wind colli- 
sion region passes closest to the WR star i.e. at i = —90°. 
As the inclination angle increases, the line of sight passes 
through less dense regions of the WR star and the opacity 
is reduced, resulting in higher observed synchrotron flux. 
For i = 0° the lines of sight to the wind-wind collision re- 
gion pass largely through the outer regions of the WR star 




Fig. 11. The 1.6-GHz intensity distribution for a model with 
D = 5 x 10 14 cm at +60° (top), +30° (middle), and -30° 
(bottom) inclination. At positive angles the OB star is closer 
to the observer than the WR star. The crosses denote the po- 
sition of the stars, with the WR star at (0,0). Each image has 
the same intensity greyscale and contour intervals (also used 
in the middle panel of Fig. 9, showing the 0° case). The op- 
tically thick part of the OB-star wind is the obvious "hole" 
(top panel), below which the synchrotron emission is seen, sur- 
rounded by the free-free emission, largely from the WR wind. 
As inclination decreases, the hole becomes less pronounced and 
the synchrotron intensity increases as the line-of-sight opacity 
through the OB wind decreases. When the WR star is closer 
to the observer (i < 0), the optically thick part of the WR-star 
wind obscures the lower part of synchrotron emission region. 



envelope and the low free-free opacity of the wind-wind 
collision region. For i > 0°, the region of the wind- wind 
collision zone occulted by the optically thick surface of the 
OB-star wind moves gradually toward the shock apex as 
i increases (cf. Fig. 11). This, and the fact that the lines 
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Fig. 12. The effect of SSA on the spectrum of the emission 
emerging from the shocked region as a function of inclination 
angle for the standard model with £ = 1CP 4 . The inclination 
angle of the axis of symmetry to the plane of the sky for the 
various spectra is shown. Dashed lines indicate i < 0° (i.e. WR 
star in front), the solid line i — 0° (i.e. quadrature), and dotted 
lines i > 0° (i.e. OB star in front). The dot-dashed line is the 
intrinsic synchrotron emission from the shocked gas, which is 
independent of inclination, plotted for reference. 



of sight to the shock apex pass through the dense inner 
regions of the OB-star wind, is the cause of the increas- 
ing absorption at low frequencies with increasing i. The 
flux when i = +90° is always greater than the flux when 
i = —90° since the optical depth unity surface of the WR 
wind is always larger than that of the OB-star wind owing 
to the higher density of the WR wind. 

The effect of inclination on the resulting intensity dis- 
tribution is demonstrated in Fig. 11, where the line-of- 
sight passes down the shock cone, and the OB-star wind 
absorbs flux from its far side. Absorption by the optically 
thick region of the OB-star wind is distinctly apparent. 
Clearly the inclination angle of a given source can have a 
large effect on the resulting flux and the spatial intensity 
distribution. 

The Razin effect only affects the generation of syn- 
chrotron photons and is independent of inclination. On the 
other hand, SSA is dependent on the path length through 
the shocked gas, and inclination angle influences its ef- 
fect on the low frequency synchrotron spectrum. This is 
shown in Fig. 12, where the intrinsic synchrotron spec- 
trum is shown. Examining the intrinsic spectrum elimi- 
nates the effect of the large changes in free-free opacity 
with inclination angle, as seen in Fig. 10. The dependence 
of SSA on the inclination angle is clearly a function of 
the path length through the shocked envelope, particu- 
larly the highest density (shock apex region) part of the 
shocked gas. Thus at high angles, this path length is small 
and the impact of SSA on the spectrum is lowest, whereas 
at quadrature the lines-of-sight through the shock apex 
region are at maximum length and the impact of SSA 
on the spectrum is highest. However, the effect is not as 
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Fig. 13. The change in thermal flux as a function of inclination 
angle using the standard model with £ = 1CP 4 . Shown are 
curves for 1.6 (solid), 5 (dotted) and 22 GHz (dot-dashed). 



dramatic a function of inclination angle as the effect of 
varying free- free opacity with inclination angle. 

When the Razin effect, SSA, and free-free opacity are 
considered together in our models it is important to bear 
in mind that both the Razin and SSA are functions of £ 
(see Sec. 3.2). As an example, when ( = 10~ 4 , the Razin 
effect is so dominant that the spectrum is influenced little 
by changes in the inclination angle. The changes that do 
occur are due mostly to the changing free-free opacity. 

In Fig. 13 we show how the thermal component varies 
with inclination angle in our standard model. The ob- 
served free-free flux is essentially constant throughout the 
range of inclinations, except when i approaches —90° and 
+90°. At these extreme inclinations, the optically thick 
parts of, respectively, the WR and OB star winds absorb 
the free-free emission arising from the bulk of the stellar 
wind of the other star and from the wind-wind collision 
region (a relatively small contribution to the total thermal 
emission (Stevens 1995)), producing a decrease in thermal 
flux. The free-free flux at i — —90° is slightly higher than 
at i = +90°, consistent with the higher density of the WR 
star wind. 



4. Modelling the radio emission from WR 147 

Having explored how the radio flux varies with separation 
and inclination in the previous section, we now turn our 
attention to the modelling of a specific system. WR 147 
is notable for being among the brightest WR stars at 
radio frequencies, and for being one of two systems in 
which the thermal and synchrotron emission are observed 
to arise from two spatially resolved regions e.g. Williams 
et al. (1997, and references therein). Furthermore, it is one 
of a handful of WR+OB binary systems where the two 
stars are spatially resolved (Williams et al. 1997; Niemela 
et al. 1998). These observations suggest a very wide sys- 
tem with a projected separation Dcos i — 0.635 ±0.020". 
At the estimated distance of ~ 0.65 kpc (Churchwell et al. 
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1992; Morris et al. 2000) this corresponds to a separation 
D <~ 415/ cos i AU. This relationship between D and i 
represents an important constraint for any models of the 
system. The inclination angle is unknown so we investigate 
models for several different values which requires different 
values of the physical separation, D, to maintain the ob- 
served angular separation. This leads to different sets of 
physical parameters, including the normalisation constant 
to fit the observed spectrum. 

4.1. Modelling the radio spectrum of WR 147 

The radio spectrum of WR 147 is perhaps the best ob- 
served of all massive binary systems, with radiomctry ex- 
tending from 353 MHz to 42.9 GHz. Reviewing the litera- 
ture, it is immediately apparent that differences in excess 
of 50% exist in the derived synchrotron flux at some fre- 
quencies (cf. Fig. 14), giving rise to uncertainty in the po- 
sition of the low frequency turnover. Skinner et al. (1999) 
estimate that at 5 GHz the synchrotron flux (S'sghz) is 
greater than the synchrotron flux at 1.4 GHz, whereas 
Setia Gunawan et al. (2001) estimate that 5i.4gh z > 
5*5GHz • These differences may be attributed to a number of 
factors that include analysis techniques, resolution effects, 
the assumed thermal contribution, and possibly tempo- 
ral variations. Though Churchwcll et al. (1992) report 
variations in radio flux approaching 50%, subsequent re- 
analysis of the same data by Setia Gunawan et al. (2001) 
suggests that any variations are much smaller (~ 15%) 
and are long-term, with little power on time scales of 
days. Skinner et al. (1999) completed their observations 
of WR 147 within a couple of days, and their repeated ob- 
servation at 1.4 GHz showed no variation. This suggests 
that, at the very least, their observed spectrum is free of 
temporal variations. In addition, it is clear there are dif- 
ficulties in determining the component fluxes when the 
sources are spatially resolved cf. Table 1 in Skinner et al. 
(1999) and Table 3 in Churchwell et al. (1992). 

These differences form the fundamental reason for con- 
flicting conclusions about the nature of the underlying 
electron energy spectrum in the current literature. While 
both Skinner et al. (1999) and Setia Gunawan et al. 
(2001) note that a free-free absorbed synchrotron power- 
law model is a poor fit to the data, Setia Gunawan et al. 
(2001) argue that the deficit in the observed flux above 
15 GHz is due to a high energy limit to the electron accel- 
eration. In contrast, Skinner et al. (1999) argue that the 
spectrum is best fit by a mono-energetic relativistic elec- 
tron distribution. At high frequencies, the spectrum from 
such a distribution falls off as P(v) oc u x / 2 e~ v . This fits 
the 15 and 22 GHz data very well, but the flux at these 
frequencies are most uncertain. Furthermore, the fit of a 
mono-energetic electron distribution spectrum to the high 
frequency fluxes depends also on the position of the low 
frequency turnover which, as already noted, is uncertain. 
How such an electron distribution results from a shock 
acceleration is unclear. 
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Fig. 14. The synchrotron spectrum of WR 147 as deduced by 
several authors - solid squares (two separate estimates from the 
same observational data by Skinner et al. (1999)), solid pen- 
tagons (Setia Gunawan et al. 2001), open circles (fluxes from 
Churchwell et al. (1992) and Contreras & Rodriguez (1999)). 
The flux at both 22 and 43 GHz is highly uncertain, with 
the latter estimated from a long extrapolation of a power- 
law fit to the lower frequency "thermal" data points to be 
« 2 mjy (Skinner et al. 1999; Setia Gunawan et al. 2001). Both 
high frequency points suggest a high frequency turn down. The 
lines are theoretical spectra from models where only free-free 
absorption in the circum-binary stellar wind envelope is atten- 
uating the synchrotron emission. Models for various inclina- 
tion angles are shown: 0, —30, —60 and —85° (corresponding 
to C = 6.92,7.23,8.08, and 13.44 x 10~ 3 respectively). The 0° 
model is denoted by the solid line. 



Using our model, we attempt to model the spec- 
trum of WR 147 assuming a power-law electron energy 
distribution. For the stellar wind parameters we adopt 
M W N8 = 2 x lO- 5 M yr- 1 , u^wns = QSOkms" 1 , 
«ooOB = 1000 km s -1 and a wind momentum ratio rj = 
0.02 (Pittard et al. 2002). The last is toward the higher 
end of reasonable estimates, and since it gives Mob = 
3.8 x 1O _7 M yr _1 favours a companion of spectral type 
around late-0 or carly-B supcrgiant. Such parameters are 
consistent with the estimated B0. 5 by Williams et al. 
(1997) based on an uncertain luminosity estimate, and the 
more recent 05-71 deduced from HST STIS spectroscopy 
(Lepine et al. 2001). The composition of the WN8 stellar 
wind is taken from Morris et al. (2000), giving X — 0.09, 
Y = 0.89, and Z = 0.016. The wind temperature for both 
stellar winds was assumed to be 10 kK, and the dominant 
ionization states were H + , He + and CN0 2+ . The ther- 
mal flux is only weakly dependent on the assumed wind 
temperature (through the Gaunt factor), if the ionization 
state is fixed. 

Values of ( are estimated for each of our models by 
matching the model to the mean of the "observed" syn- 
chrotron 4.86 GHz fluxes, which we take to be 14.1 ± 
0.3 mjy. We then determine the spectrum between 0.2 
and 80 GHz assuming this constant value of £. As the in- 
clination of the model increases, the binary separation D 
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Fig. 15. Model synchrotron spectra when SSA, the Razin ef- 
fect and free-free absorption in the circum-binary stellar wind 
envelope are included, for various inclination angles: (solid), 
±30 (dotted), and ±60 (dashed), with corresponding values of 
C = (7.03, 7.32, 8.12) x 10~ 3 . The data points are those shown 
in Fig. 14. 

has to increase in order to satisfy the observed constraint 
on Dcosi, forcing the need to increase £ to maintain the 
5 to 15 GHz flux level. 

To fit the thermal flux with our chosen value of M 
we require the winds to be clumped. This has some 
precedent: Morris et al. (2000) note that the observed 
thermal flux implies that M — 7.5 x 10~ 5 M yr _1 if 
the wind is assumed non-clumpy, but find that models 
with a homogeneous wind cannot match the 1R helium- 
line profiles observed with ISO. Their derived volume- 
filling factor of / ~ 0.1 implies an actual mass- loss rate 
M ~ 2.4 x 10~ 5 M Q yr _1 . We have therefore multiplied 
£^ and from Eqns. 1 and 2 by a factor 1// to simulate 
clumping in the cool stellar winds. We assume that both 
the WR and O-star winds are clumped by the same factor. 
Also, it is assumed that the clumps are destroyed as they 
pass through the shocks into the wind collision zone, and 
therefore do not increase the free-free emission from this 
region. We adopt / = 0.134 2 to ensure that the thermal 
flux at 22GHz remains less than the total flux. 

Fig. 14 demonstrates that the low-frequency turnover 
could potentially be accounted for in our models by free- 
free opacity of the circum-binary stellar wind at interme- 
diate inclination angles, with the line of sight into the 
system through the WR wind 3 , but not at 0° or incli- 
nations angles for which the lines-of-sight to the shock 

2 Interestingly, this implies that the mass-loss rate of the 
equivalent homogeneous wind is M — M a C t ua i / V7 = 5.5 x 
10~ 5 M Q yr _1 . While this is below the observationally deduced 
value using a single star model, it is consistent with the 13- 
38% enhancement in S v expected from a binary with a wind 
momentum ratio n = 0.1 (Stevens 1995). 

3 Dougherty et al. (2003) erroneously report a fit to the low- 
frequency spectrum of WR 147 using only free-free absorption. 
Their fit was attained using SSA, in addition to free-free ab- 
sorption. 



0.1 1 10 100 

Frequency (GHz) 

Fig. 16. The spectrum of WR 147. The synchrotron data are 
those shown in figs. 14 and 15 and represented by solid data 
points. The total flux as deduced by the same authors are repre- 
sented by the hollow data points. The lines are the total (solid) , 
thermal (dotted) and synchrotron (dot-dashed) spectra of the 
i = 0° models with £ = 7.03 x 10~ 3 , where SSA, Razin and 
free-free absorption are included in the radiative transfer cal- 
culations. The volume-filling factor / = 0.134. Models at ±30° 
and —30° are also consistent with the data. 

apex pass through the densest parts of the WR wind ie. 
—85°. The models do not fit the 1.4 GHz fluxes particu- 
larly well. Setia Gunawan et al. (2001) arrived at a dif- 
ferent conclusion: that free-free absorption alone was suf- 
ficient to account for the turnover, based on results from 
a synchrotron point-source model. We attribute our dif- 
ferent conclusion to the extended nature of the region of 
shocked material in our model; although the line-of-sight 
to the shock apex may be optically thick, lines-of-sight to 
the outer regions of the wind-wind-collision region may re- 
main optically thin. Maximum free-free absorption of the 
low frequency synchrotron emission occurs around incli- 
nation angles of —60°. For still lower inclination angles, 
the observed flux is higher despite increased attenuation 
to the shock apex. This is because increased £ leads to 
synchrotron emission further from the shock apex. The 
only manner in which free-free opacity alone could pro- 
duce a sufficiently strong low-frequency turnover would 
be if the size of the synchrotron emission region were re- 
duced, such that the greater attenuation near the shock 
apex became more significant to the observable emission. 
Such a reduction in the size of the emission region could be 
mimicked in an ad-hoc fashion in our model by evolving £ 
away from the shock apex. Certainly, as the shocks become 
more oblique moving away from the shock apex, the parti- 
cle injection and acceleration efficiency is lowered (Ellison 
et al. 1995, 1996). The energy spectrum is also evolving 
since IC cooling is ongoing as the electrons advect away, 
though this predominantly effects the higher energy elec- 
trons, and coulombic cooling may be more important for 
the lower energy electrons. 

Models that include the influence of both SSA and the 
Razin effect in addition to free-free absorption are more 
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Table 1. Values of density, equipartition field, and peak tem- 
perature near the shock apex for the models shown in Fig. 15. 



Inclination 


0° 


30° 


60° 


D x lO" 15 cm 


6.21 


7.17 


12.4 


n e (WR) x 1(T 4 cm- 3 


1.8 


1.35 


0.45 


n e (OB) x 1(T 4 cm" 3 
B 1 mG 


3.3 


2.5 


0.83 


4.1 


3.6 


2.2 


C x 10 3 


7.03 


7.32 


8.12 


T max (WR) K 




2.4 x 10' 




T max (OB) K 




1.4 x 10 7 





The B-field at the shock apex is approximately the same in 
both the shocked WR and OB regions. 

successful in fitting the low-frequency data. In Sec. 3, both 
these effects were shown to be significant, depending on 
the value of £, and it seems unrealistic to exclude such fun- 
damental processes from our models. The resulting spec- 
tra are shown in Fig. 15. As before, values of £ have been 
chosen for the different models based on matching the 5- 
GHz data point, and are ~ 1%, similar to the values used 
in models of SNR (cf. Mioduszewski et al. (2001)) and 
of an order of magnitude expected for shock-accelerated 
electrons (see Sec. 2.2). As in the case of free-free opac- 
ity alone, as the separation increases, ( increases and the 
amount of emission arising far from the shock apex in- 
creases. Thus the low-frequency turnover from optically 
thin to optically thick emission moves to lower frequen- 
cies with increasing separation. The difference between 
the spectra for positive and negative inclinations of the 
same absolute value is due to the different free-free opac- 
ity through the OB-star wind as opposed to the WR-star 
wind. The values of density, equipartition field, and peak 
temperature near the shock apex for the models shown in 
Fig. 15 are given in table 1. 

In Fig. 16 the resulting total, synchrotron and thermal 
spectra from one of our models (i — 0, and including SSA, 
the Razin effect and free-free absorption) is shown against 
the observed spectra. Though this is not the best fit to the 
data in a formal sense, our model fits the data reasonably 
well given the approximations and assumptions which it 
contains. In the model shown in Fig. 16, the total flux 
in the 5 to 8 GHz range is underestimated, and the syn- 
chrotron emission does not turn down somewhere around 
10 to 20 GHz to account for the observed high frequency 
data. We could account for the shortfall in the 5-8 GHz to- 
tal flux by increasing the thermal emission by a few mJy. 
This could be achieved by simply increasing M /voo in the 
stellar winds. However, this would require that the syn- 
chrotron emission at frequencies higher than <~ 10 GHz 
be lower than in the current models, which may well be 
the case if IC cooling of the highest energy electrons is 
taken into account, as explained below. 

Fig. 16 lends considerable support to the suggestion 
that WR stars with continuum spectra of spectral index 
less than +0.6 are synchrotron emitters e.g. Dougherty 
& Williams (2000). The total flux spectrum as shown 
in Fig. 16 could be fit with a power-law spectrum of 



spectral index of <~ +0.3, which could be interpreted 
as arising from a partially optically thick stellar wind. 
However, it is well-established that the stellar winds of 
WR stars have radio continuum spectra with spectral in- 
dices +0.6 -> +0.8 e.g. Williams et al. (1990b); Lcitherer 

6 Robert (1991), and spectral indices less than <~ +0.6 
are not expected from stellar winds alone. This suggests 
that spectral indices less than +0.6 are a strong indica- 
tor of the presence of a synchrotron component. This is 
especially useful in those cases e.g. WR86, where the syn- 
chrotron and thermal emission regions are not resolved as 
separate components, yet a continuum spectrum with a 
low spectral index can be observed. 

We are unwilling to draw much of a conclusion related 
to inclination angle from these models of the spectrum 
of WR 147. For models with wider separations i.e. higher 
inclinations, the 1.4 GHz data is overestimated substan- 
tially, which hints at an inclination angle somewhere be- 
tween 0° and 30—40°. The high end of this estimate is con- 
sistent with the poorly constrained inclination estimates 
from Williams et al. (1997) and Contreras & Rodriguez 
(1999). However, any conclusion related to inclination an- 
gle drawn from fitting the spectrum has to be tempered 
due to the heavy reliance on the accuracy of the 1.4 GHz 
data points, and the poor precision of the 352 MHz data 
point. 

Another striking feature of Figs. 14, 15 and 16 is our in- 
ability to generate a high-frequency turnover somewhere 
around 10 to 20 GHz, as suggested by the data points. 
This is a consequence of using Eqn. 6 to describe the syn- 
chrotron power, where it is assumed that 7 max is infinite. 
IC cooling first determines the maximum energy that can 
be achieved by the acceleration processes. As the relativis- 
tic electrons are advected away from the acceleration site 
by the flow, the highest energy particles suffer crippling 
losses from IC cooling as this is no longer balanced by 
mechanisms supplying rapid energy gain. Hence 7 max will 
have a finite value. For WR 147, the rate of energy gain by 
lst-order Fermi acceleration is equal to the rate of energy 
loss by IC cooling when 7 = 7 max ~ 10 6 . From Eqn. 16, 
electrons with 7 = 10 6 are cooled to 7fl ow ~ 400 within 
a time ifl ow when r w tob and L„ = 10 5 L Q . The high 
frequency turn down seen in the observational data oc- 
curs at <~ 10 — 20 GHz, which for B <~ 3 mG suggests 

7 ~ 1200, which is within a factor of a few of the value es- 
timated above. Thus, we conclude that the high frequency 
turn down is an expected consequence of IC cooling. This 
break in the energy spectrum of the relativistic electrons 
will produce a corresponding break in the spectrum of 
the IC photons. For relativistic particles, the final energy 
(Ef) of an IC photon is related to its initial energy (Ei) by 
Ef ="f 2 Ei. For an OB-type star, the stellar photon distri- 
bution peak occurs ~ 10 eV. This implies that the IC pho- 
tons resulting from scattering off electrons with 7 ~ 10 3 
or greater will have energies in excess of 10 MeV - the hard 
X-ray and gamma ray regimes. Such high energy photon 
production in CWBs has been examined by other authors 
e.g. Benaglia & Romero (2003, and references therein). 
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An additional point is that since dj oc 7 2 , the higher en- 
ergy electrons cool most rapidly. This causes the electron 
distribution to "bunch-up" near 7fl OW) an d may help to ex- 
plain why a mono-energetic fit works well. A model with 
an evolving value of 7 will be explored by Pittard et al. (in 
preparation) . 

4.2. Modelling the MERLIN images of WR 147 

The spatial distributions of radio emission at both 5 and 
1.6 GHz from the model in Fig. 16 are shown in the top 
panels of Fig. 17. At 5 GHz the peak intensity of thermal 
emission from the WR stellar wind and the synchrotron 
emission from the collision region are similar, and the OB- 
star companion is clearly visible. At 1.6 GHz the syn- 
chrotron emission is much brighter than the WR star stel- 
lar wind, and the OB-star wind is barely visible. 

To appreciate how these model intensity distribu- 
tions compare with the MERLIN observations shown 
in Williams et al. (1997), we used the AIPS subroutine 
UVCON to generate visibilities appropriate for a MERLIN 
"observation" of our models. In addition to the array co- 
ordinates, system noise estimates were included in the 
calculation by including the performance characteristics 
of the MERLIN telescopes e.g. antenna efficiencies, sys- 
tem temperatures etc. The resulting visibilities were then 
imaged and deconvolved using the same procedure as 
in Williams et al. (1997), giving the "simulated" observa- 
tions shown in the lower panels of Fig. 17. The remarkable 
similarity between these images and the observations pre- 
sented in Williams et al. (1997) is striking. However, here 
the 5-GHz peak intensity of the collision region is a little 
lower than that of the WR stellar wind, and the WR star 
is a little too bright at 1.6 GHz. The WR star also appears 
to be a little more compact than shown in Williams et al. 
(1997), but we attribute these minor differences to our 
simple spherical, isothermal model of the stellar winds. 
Note that the OB star is not visible in the simulated ob- 
servations. This is because the observations are essentially 
the intensity distribution (the upper panels) convolved 
with the interferometer beam and the emission from the 
OB star is then sufficiently dispersed that it is no longer 
visible. 

One of our hopes in examining the simulated images 
was to help constrain the inclination angle. The observa- 
tions shown in Fig. 17 are for an inclination of 0°. At an 
inclination angle of 30°, the synchrotron emission is spread 
out over a much larger area and has a much lower surface 
brightness than in the observations for i = 0°, and shown 
in Williams et al. (1997). As from the spectral modelling, 
this leads us to believe that the inclination of the system 
is quite low, between 0° and 30°. This is contrary to the 
conclusions of Pittard et al. (2002) where larger inclina- 
tion angles were required to recover the extended X-ray 
emission. However, an alternative explanation for this is 
that the stellar X-ray emission is extended on larger scales 
than previously thought (see Skinner et al. 2002). 



5. Summary and Future Directions 

In this first paper we have modelled the thermal and syn- 
chrotron emission from early-type binary systems with 
a strong wind-wind collision. We have used appropriate 
simplifying assumptions to make this initial investiga- 
tion tractable. In particular, we have generated models 
of wide systems without considering IC cooling, and have 
assumed that the relativistic particle energy distribution 
is a power-law up to infinite energies and is spatially in- 
variant throughout the shocked gas. We have also assumed 
that the magnetic field is highly tangled so that the syn- 
chrotron emission is isotropic. 

We have demonstrated the importance of considering 
extended emission and absorption regions, as opposed to 
treating these in a point-like manner as in previous work. 
We have also demonstrated the importance of synchrotron 
self-absorption, the Razin effect, and the free-free opac- 
ity of the shocked gas in these systems. With our rather 
simple model we have been able to model both the spec- 
trum and the spatial distribution of radio emission from 
WR 147 remarkably well. However, our simple model fails 
to address the apparent high frequency cut-off in the syn- 
chrotron spectrum, due to the current lack of a high energy 
cutoff in the relativistic electron energy distribution. We 
note that this cutoff, which suggests a characteristic value 
of 7 ~ 1200, can be naturally explained by IC cooling. 

Clearly if we are to realistically model the observed 
radio emission from CWBs we will need to calculate the 
evolution of the relativistic particle energy spectrum as 
the particles are advected downstream. This scenario will 
be addressed in a follow-up paper (Pittard et al. , in prepa- 
ration) where we plan to introduce a finite maximum en- 
ergy for the relativistic particles and allow this to change 
along the shock front, as expected for variable IC cooling 
by photons from the OB star. The inclusion of IC cooling 
will enable the calculation of the evolution of the energy 
spectrum of the relativistic particles downstream of the 
shock acceleration region. Similarly the effect of Coulomb 
cooling on the lower energy relativistic particles will be 
incorporated. The emission spectrum of mono-energetic 
particles can then be convolved with the energy spectrum 
of the particles to determine the total emission. 

With these additional mechanisms it will be possible 
to investigate the high frequency cut-off, such as that ap- 
parent in WR 147, and to more accurately model the low 
frequency cut-off when coulombic cooling is important. 
Addressing these processes is vital if we are to model 
smaller binary systems where the electron energy spec- 
trum is highly variable, largely due to the IC cooling 
from stellar UV photons. This is certainly the case in the 
proto-typical CWB system WR 140, where it is clear the 
relative timescales for shock acceleration of electrons, IC 
cooling, coulombic cooling and adiabatic expansion are 
varying dramatically throughout the 7.9-yr orbit. At the 
very least, these effects, along with radiative shock cooling, 
need to be considered in any realistic model of WR 140. 
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Fig. 17. Intensity distributions at 4.8 (left) and 1.6 GHz (right) of our WR 147 model at 0° inclination. The top panels are the 
model intensity distributions at each frequency. The greyscale and contour range are the same in both images. The stellar wind 
of the OB-star companion is visible, particularly at 5 GHz, along with the WR star wind and the wind-wind collision region. 
The large contour range was used to show that the stellar wind of both stars extends far beyond the greyscale range. The lower 
images are simulated MERLIN observations generated using the model intensity distributions shown at the top, and the same 
u — v distribution and beam sizes as Williams et al. (1997). The beams are circular, of diameter 57 and 175 mas at 5 and 1.6 
GHz respectively. The contour levels at the two frequencies are also the same as those in Williams et al. (1997). The similarities 
to the observations of Williams et al. (1997) is striking, most particularly the emission from the wind- collision region. 
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Appendix A: Radiative Transfer 

The equation of radiative transfer in a medium in which 
emitting and absorbing material is present, in the absence 
of scattering, can be written 

= £„ — a v I v (A.l) 

as 

where e v and a v are the emission and absorption coeffi- 
cients at frequency v, s is the co-ordinate along a line of 
sight, and I v is the intensity of the radiation field. This 
has a solution, at a given frequency v 

I(s) = I(0)e- T(s) + f e(s')e- r(s '^ s) ds' (A.2) 
Jo 

where 

T (s) = f a(s')ds' and r(s' -» s) = f a{s")ds" . (A.3) 

JO Js' 

Consider a line of sight along which a and e are piecewise 
constant, that is 

e = Ei and a — on for (i — l)As < s < iAs (A. 4) 

where As is a finite increment of path length and i = 
1 . . . n. Solving the integrals in the solution to the radiative 
transfer equation (Eqn. A.2) we can write the emergent 
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intensity from an increment of path length (i — l)As < 
s < iAs as follows 



Ii-ie 



-a.As 



iAs 



(i-l)As 



(A.5) 



/ e l e-^ s 'ds' (A.6) 

J(i-l)As 

(A.7) 



h = k- ie - a > As + - (1 - e-^ As ) . 

Starting with 1\ , and applying this scheme to each incre- 
ment of path length in turn, we can derive the emergent 
intensity /„ along the line of sight. In general, the incre- 
ment of path length As need not be constant, as long 
as the emission and absorption coefficients are constant 
within a given As. 

This scheme forms the basis of the code. Images are de- 
rived by calculating the emergent intensity on a number of 
lines of sight which pass through a grid of emission and ab- 
sorption coefficients defined in axis-symmetric cylindrical 
polar coordinates. The geometry of the problem is shown 
in Fig. A.l. It is assumed, without any loss of generality, 
that the line of sight which passes through the origin of 
the coordinates on the model grid also passes through the 
origin of the coordinates on the plane of the sky, and that 
the projection of the line of sight onto planes of constant 
z is parallel to the line cf) = 0. 

The cell j, k on the model grid is bounded by 

(j - l)Ar < r < jAr and (k - l)Az < z < kAz (A. 8) 

where Ar and Az are, respectively, the width and height 
of a cell in the r, z plane. The grid consists of N r cells in 
the r direction by N z cells in the z direction. The domain 
covered by the grid is therefore 



r < i?;0< z < Z 



(A.9) 



where R — N r Ar and Z — N z Az. The z axis of this grid 
is inclined at an angle 9 to the plane of the sky, on which 
we define a co-ordinate system u, v. A line of sight through 
the grid is therefore defined by u, v and 9. All lines of sight 
are normal to the u, v plane. 

It is convenient to define a point along a line of sight 
in terms of its azimuthal angle <p in the frame of the grid 
of emission and absorption coefficients. The co-ordinates 
on this grid of the point (j) along the line of sight u,v,6 
are given by 



r = 



sm< 

v 



z = 



3S6» 



tan 9 
tan 4> 



(A.10) 
(All) 



The point, </>o, at which the line of sight enters the grid is 
given by 



»o = mm 



injsfiT 1 ^tan -1 (-^ sin 6»)j . (A.12) 
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Fig. A.l. Geometry for solving a line of sight through a grid 
in cylindrical polar coordinates 

The value of 0o can be used to determine the indices (j, k) 
of the first grid cell which encounters the line of sight, 
using Eqn. A. 8, A. 10, and A. 11. The indices of the next 
cell are determined from the minimum of 4> r and <j> z , which 
are the values of cf> corresponding to the boundaries with 
the next cells in the r and z directions respectively. The 
path length through the current cell is then given by 



As = uy/l + tan 2 9 ( — 1 — 
\ tan 0i 



1 



tan< 



(A.13) 



where 4>\ is the point at which the line of sight crosses the 
next cell boundary. 

For a given line of sight, we now have the grid indices 
of the first cell which encounters the line, the indices of the 
next cell, and the path length in the current cell. Having 
obtained the emission and absorption coefficients from the 
appropriate grid cell, we can update the intensity on the 
given line of sight using Eqn. A.7. This process is repeated 
for a line of sight until it leaves the grid (i.e. if i > N r or 
j > N z ). 



